home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / fractal / kaos.lha / fixptlib / trkman_ode.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-01-14  |  3.1 KB  |  116 lines
  1. /*
  2. ----------------------------------------------------------------------
  3.      Compute a set of manifolds for saddle equilibria of ODE
  4. ----------------------------------------------------------------------
  5.      compute manifolds of ANY given set of hyperbolic fixed points
  6.      with eigenvalues and eigenvectors
  7.      iskip : number of skips between each data
  8. Input:    sevec[n][n]: eigenvectors for n eigenvalues
  9.     seval[n]: eigenvalues of an equilibrium
  10.     xe[n]: coords of an equilibrium
  11.     mu: # of maximum steps along unstable manifold
  12.     ms: # of maximum steps along stable manifold
  13.     muf : # of finer division of the unstable manifold
  14.     msf : # of finer division of the stable manifold
  15.     iskip: # of skips before drawing
  16.     delman: initial distance from an equilibrium 
  17.     n: phase space dimension
  18.     r: unused, always set to 0 within the subroutine (kept to retain the
  19.     same syntax as trkman_map.c)
  20. ----------------------------------------------------------------------
  21. NOTE:    -Only works for r=0; (equilibrium point of ODE)
  22.     -Internal variable: man_index: -1=complex, 0=unstable, 1=stable
  23. ----------------------------------------------------------------------
  24. */
  25. #include "../include/x11r2_kaos_def.h"
  26.  
  27. #define RK4 1    /* Runge-Kutta integrator */
  28. #define DELFRAC 0.1
  29. void trkman_ode(sevec,seval,xe,r,ms,mu,msf,muf,iskip,delman,n)
  30. int r,ms,mu,msf,muf,iskip,n;
  31. double **sevec,**seval,*xe,delman;
  32. {
  33.     int i,ic,jc,mc,mend,man_index,icnt,color;
  34.     double fabs(),t_step,delx,*x,*xp,*x1,*x2,*dvector();
  35.     extern int stop,int_algorithm,dot_size;
  36.     extern double time,time_step;
  37.     extern char string[];
  38.  
  39.     r=0; /* r should be 0 */
  40.     x1 = dvector(0,n-1);
  41.     x2 = dvector(0,n-1);
  42.     x = dvector(0,n-1);
  43.     xp = dvector(0,n-1);
  44.  
  45.     for(ic=0;ic<n;ic++){    
  46.         if(seval[ic][1]!=0) /* if complex eigenvalue quit*/
  47.             man_index = -1;
  48.         else {
  49.             if(seval[ic][0]>0){
  50.                 man_index = 0;
  51.                 color = Blue;
  52.             }
  53.             else {
  54.                 man_index = 1;
  55.                 color = Red;
  56.             }
  57.         }
  58.         if(man_index== 0){
  59.             t_step = time_step/muf; /* temporary time step */
  60.             mend = mu * muf;
  61.         }
  62.         else if(man_index == 1){
  63.             t_step = - time_step/msf;
  64.             mend = ms * msf;
  65.         }
  66.  
  67.         sprintf(string,"Computing manifold for ev[%d]...",ic);
  68.         printf("%s\n",string);
  69.         printf("(Ev=%g %g)\n",seval[ic][0],seval[ic][1]);
  70.         printf("(Evec %d =%g %g)\n",ic,sevec[ic][0],sevec[ic][1]);
  71.  
  72.         for(jc=0;jc<2;jc++){
  73.             if(man_index==0){
  74.                 printf("Unstable manfold for ev[%d]...\n",ic);
  75.             }
  76.             else if(man_index==1){
  77.                 printf("Stable manfold for ev[%d]...\n",ic);
  78.             }
  79.             if(jc==0){
  80.                 for(i=0;i<n;i++) x[i] = xe[i]+delman*sevec[ic][i];
  81.             }
  82.             else {
  83.                 for(i=0;i<n;i++) x[i] = xe[i]-delman*sevec[ic][i];
  84.             }
  85.             icnt=0;
  86.             for(mc=0;mc<mend;mc++){
  87.                 for(i=0;i<n;i++){
  88.                     delx=(x[i]-xp[i])/seval[ic][0];
  89.                     x2[i]=x[i]+delx;
  90.                     x1[i]=x1[i]- DELFRAC*delx;
  91.                 }
  92.                 /*use Runge-Kutta 4th order */
  93.  
  94.                 int_onestep(x1,x2,x,&time,t_step,n,RK4);
  95.                 for(i=0;i<n;i++){
  96.                     xp[i] = x[i];
  97.                     x[i] = x1[i];
  98.                 }
  99.  
  100.                 if((icnt % iskip)==0){
  101.                     (void) draw_record_pwf(x1,color,-1,dot_size,0,1);
  102.                 }
  103.                 icnt++;
  104.                 if(stop){
  105.                     goto done;
  106.                 }
  107.             }
  108.         }
  109.     }
  110. done:
  111.     (void) free_dvector(x1,0,n-1);
  112.     (void) free_dvector(x2,0,n-1);
  113.     (void) free_dvector(x,0,n-1);
  114.     (void) free_dvector(xp,0,n-1);
  115. }
  116.